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ABSTRACT 


The velocity field and pressure distribution between stationary 
and rotating finite disks are obtained by a finite difference solution 
of the Navier-Stokes equations and the continuity equation using an 
explicit scheme developed by Chorin [1]. These numerical results 
are compared with the results of "infinite disk’ solutions. In the 
range of interest for viscometer applications, the velocity and pres- 


sure are adequately represented by the infinite disk solutions. 
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Symbol 


P 


P 
(r, 6, 2) 
(R, 6, 2) 


R 


R 


TABLE OF SYMBOLS 


Definition 
Artificial density 
Number of iterations 
Pressure 
Dimensionless pressure, D/6 
Cylindrical coordinates 
Dimensionless cylindrical coordinates 
Maximum radial distance 
Reynolds number 
Time 
Dimensionless time 
Radial Valocity 
Dimensionless radial velocity 


Dimensionless radial velocity from 


the infinite disk solution 
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Tangential velocity 


Nimensionless:tangential velocity 
Axial velocity 

Dimensionless axia) velocity 
Maximum axial distance ” 
Artificial compressibility factor 
Ratio, Z/R 

Viscosity 


Kinematic viscdsity, p/p 


Density vital asker aah 


Angular velocity * 


T* 


I, INTRODUCTION 


The problem of practical interest to be solved can be stated as 
follows: An incompressible Newtonian liquid is contained between an 
upper, stationary disk and a lower, rotating disk. Determine the 


velocity field and pressure distribution of the fluid in motion. 


— Stationary Disk 
—Liquid-Air Interface 
—— Rotating Disk 


2 


Figure 1. Physical Configuration 


By presenting accurate ntitpemteae solutions to the Navier-Stokes 
equations, this paper makes possible an effective determination of the 
nature of inertial secondary motions in a rotating disk viscometer 
having disks of finite radii. Previous considerations of the problem 
have assumed that the disks are of infinite radii, resulting in solutions 
involving numerical and analytical approximations [2, 3]. A compari- 
son will be made of the velocity distributions obtained with disks of 


finite radii and those having radii assumed to be infinite. 


F 


Il, THE EQUATIONS OF MOTION 


The equations governing flow complet of those ot Nae wakes 
together with the equation for continuity. “These, as s used in this 
paper, are for a y ici oendingte system fe 6, z) al, rota- 

t 


tional symmetry and incompreseibility saanened: The saidetious are 


written as follows [4]: : 


ele aE aa 


Here, u, v, andw xeyeencee respectively, velocities in the 
on) ae 2 {4 Je 15 5 ; 
radial, ba ag th and axial directions. p is pecs P is press sure; 
USO Be a ¢ 
and py is Kinematic eter, ee ‘The 5 seer ig of wei disk is R; the 


ee 7 =3F3 rf ee 


‘spacing between them is Z and the angular velocity of the Potating 


; 
ifs 


disk is 9. 
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‘The above equations are normalized-according to the following 


procedure: 
Primary Variable Normalizing Quantity Normalized Variable 
u, Vv, W RQ U, V, W 
r | R- R 
z Zz Z 
t Z/RQ T 
P — -p R99? P 


The normalized equations are 


au 2 

ort ACo ee V)+w id as: 
af ye fu ,1 ww _ oy, 20) 
a. aR? + R aR Ret az 


aw aw aw aP. i afew 13 ew - 
—— —, ee -—— + 
at tA VaR t WS az ¢R [a Set tes )*az] 


1] 


(5) 


(6) 


(7) 


(8) 


The boundary conditions used are shown in Figure 2. 


< 
it 
=) 


IF 
© 


oY 
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" 
° 


U-W=0; V=R 


Figure 2, Boundary Conditions 


- 


a 


The Fotating plate is at Z = 0 and the Hanioantye one at Z = 1. 
\ poe) 


R = 0 is the axis ‘of rotation and R= is the liquid-air interface which 


is assumed to be planar. 


(y} Hey: : ES : a ee vu a F , 
~~ The Reynolds Number, Ro represents the méasure of inertial 


stresses one to viscous stresses. Viscometers are usually 
4 anaiysed by paid ds Bee. the “inertialess" limit ‘of R. +0, in which 
case the flow consists of circular wicca s3r Viegorieter’ gen- 
2 erally operate with a gap-to-radius ratio, 2, of less than 0; 05 and 


with the product A R, less than unity. 
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Il. NUMERICAL PROCEDURE 


To solve the Navier-Stokes equations, Chorin's explicit tech- 
nique is used. An artificial compressibility factor, 6, is introduced 
into the equations of motion. Hence the continuity equation is re- 


placed by 


oD - _ UD. 3av\_ ow 
aT a (gids 3Z (9) 


where P = D/6 is an artificial equation of state. 


Equations (5), (6), (7), and (9) are put into finite difference form 
by means of a Dufort-Frankel scheme, (details are obtained from 
Chorin [1]), and the resultant equations appearing below are then 


utilized in the computing procedure: 
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U (TT) = AAT y Cucrt) - U (RO )] 


{ -4a2 


+2AAT = 


- 2S w Cu(2t) - u (z- )] 
~ AAT {pirt) > D(R-)] 
CORAL fous sii . 
aes et ae 


ROR R 


Ti 


- 24aT u 
R, ® 


2 
+ Raz Lu (z+) + U(z-) - U(T-)) 


eur} / {1+ PSF + a az} i 


AR? * 


_ AAT 
AR 


V (Tt) = { uU LV (Rt) - V(R7)] 


- 2’\AT = 


-& w Ev(zt) - v(z-)] 


2 
+R aR [v (Rt) + V(R7) -V(T")] 


+i 4T CV iR+) - v(R-)] 
; R 


R_ AR 
e 


+R Ae Lv (z+) + V (2-7) - V(T-) 1 


2X7 AT so. 2AT 
az? f 


| + V(t) } / {) * RAR? *R AZ (11) 
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W (Tt) = 


; _ AAT 
os “AR U [Ww (Rt) - W(R-)) > 


AT : 
-az WV { w (Z+) - W (Z-)] 


3 Be [p (zt) + D(z-)] 


i rN AT 
+ Rane Cw (Rt) +.W (R-) - W(T-)) 


42or CwiRt) - wiR-)] 
RAR R= 


2AT 
* Raz Cw (z*) + W(2-) - W(T")) 


wor }/ (us BBE wy Raz) (12) 
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U 
«is ¥ 
D (TT) 2X AT R 


_NAT eee 
aR LU(RT) - U(R)] 


AT % 2 2 
ie [w (Zt) - W(z7)] + D(T-) (13) 


where U (Rt) = U(R+AR, Z, T); 


U (R~) = U(R-AR, Z, T); ete. 


Equations (10) through (13) are applicable only at interior 
points. For the various boundaries, the following equations were 


developed using for the purpose one-sided differences. 


At the centerline, (R = 0), 


U-=0 (14) 


V=0 (15) 


+). j AT +) - ~ 
WARS) a4 az WCW (2*) W (Z7)) 


-~ At 
6 AZ 


[D(zt) - D(Z-}j] 


202 aT 2 
+R aR? [2 Ww (Rt) - W(T-)) 
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2AT 


+ RAZ? [Ww (z+) + W (2) - WwW (T~)] 


¢ 7 AT 
sworn }/ {1+ Baer + gaz | (16) 


U (RT) 


~ SE Cw (zt) - w(z-)) + DIT). (17) 


At the liquid-air interface, (R = 1), 


U aiOasinitit bebie-one | _ 08) 


vert) = { - SS ow ivizt) - v2) 


EAs [2 Vv (R--) - 4(R7) + V(T")] 


2AT = 3 
+ Raz Cv (zt) + V(Z7) - V(T")) 


how 4R 


ee, Cae nae poor} (19) 
4 ts wr e 


3) 


w (Tt) = { - 85 w [w (zt) + w(z7)] 
_ AT +) 
o (pizt) - piZ-)] 


2 
+ > [2 WwW (R--) - 4W(R-) + W(T-)] 


+ AZ Lw (zt) + W(Z-) - W(T-)] 


a ar 2 2 ' 
two} /{1-RBS 28} eo 
e e 
2AAT 
+ = — 
D (TT) AR U (R7-) 
_AtT +) _ - ~ 
AZ [w (Zt) - W(Z7)] + D(T). (21) 
Along the rotating disk, (Z = 0), 
U =0 : (22) 
V=R (23) 
W= 0 (24) 


D(T+) = (= 3 3S . Lwlzt) 3) WZ) + D(T>), (25) 


Along the stationary disk (Z = 1), 


L 


U=0 4 (26) 

V=0 (27) 

Ww-=0 . (28) 

pitty = - 24% (w(z) - w(z-)] + Dit). (29) 
)¥ ‘ / : . 


3 
The above equations(to obtain steady state solutions) were pro- 


grammed for a ie 7% computer, Such solutions are considered to 


Cae) 


have been Gellteved heh D(Tt) and p(t-). are approximately equal. 


A copy of the computer program, coded in Fortran IV, is contained 


in Appendix A. oi Teh es 
St) oo q 1 
tA i 

ij A 
( 7 } YW (*si 8 a 
vx 

(Q v) itsiot tf 

{ na P Q 

{cS 1 ; 

Ww 
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IV. STABILITY 


If the flow under consideration is to remain incompressible 


and therefore stable, Chorin asserts that 6 must be chosen so that 


ot R, Max (U? + V? + w?)2 <1 
R, Z. 


¢ 


This is equivalent to keeping the artificial Mach Number less than 
unity. 
To insure that motion can be adequately described in the con~ 


fines of a given-grid scheme AT must be chosen so that 


AT < .437 Min (AR or AZ) ot 
R,Z | y: 


Numerical results show that the above conditions are not sufficient 


for stability in some cases and too restrictive in others. 
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V. NUMERICAL RESULTS AND DISCUSSION 


iA: i 
The secondary flows found in the R-Z plane are shown for 


A= 1.0, .50, .20 and. 10 in Figures 3 through 6.. (A R..= 1.0 in each 


case.) Computer progran: outputs for velocities and pressures cor- 


responding to these cases are given in Appendix B. 


- Once an adequate 6 and AT are chosen, ‘convergence is achieved 


t 


in 2000 to 4000 iterations. This requires 20 to 40 minutes of com- 


- . puter time. The following table contains 6 and AT values for which 


' divergence occurs: 


e 
.75 1.0 
75 1.0 
£75, 1.0 
1.0 1.0 


rbd ts 7 
Chorin's 


Trial Values... Necessary Conditions 
6 AT i AT | 
, N sss! 1 Fs 
- 04 001. 5 1,0 - 022 
. 09 é 001 he 0 e 022 
eat { ont 4 Fa f 

- 16 -001 1,0 - 022 

s ba k 
- 04 - 0005 1,0 - 022 
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Figure 3. 


Secondary Velocity Field 
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Listed below are 6 and AT valucs for. which there are convergent 


solutions: 
Chorin's 
Trial Values Necessary Conditions 

6 R, 5 AT 6 AT 
1.0 1.0 .04 00025 1.0 022 
75 1.0 .04 .00075 1.0 . 022 
£75 1.0 .04  _.00075 1.0 . 022 
75 1.0 £09. 00075 1.0 022 
75 1.0 .16 0005 1.0 . 022 
75 1.0 .16 .00075 1.0 . 022 
.75 —-:10.0 .04 0005 01 . 002 * 
50 1.0 09. 001 10 022 
. 50 2.0 .04 =. 001 625 011 
. 20 1.0 .09-. 002 1.0 . 022 
. 20 5.0 .04  .0005 . 04 . 004 
. 20 5.0 .04 .001 04 . 004 
.10 1.0 .04 =. 001 1.0 . 022 

10 = 10.0 .04 =. 001 01 002 * 


- 10 5.0 . 04 001 . 04 . 004 


*The necessary conditions are violated, yet a convergent solution 
is obtained. 
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At least 20 grid spaces are necessary in R and Z for reasonable 
results. Numerical oscillations are produced in the Z direction when 


only ten grid spaces. are used. 


For the small secondary flows, true convergence is very 
difficult to discern with this numerical technique. As an example, 
forr =.1 and R. = 5.0 convergence seems to have wean achieved in 
2000 iterations. But when the results are plotted, numerical oscil- 
lations in U are evident. Nowavee: after 4000 iterations these oscil- 


lations have vanished and the results appear to be accurate when 


compared with the infinite disk solution. 


Rate of convergence is independent of starting values unless 
accurate values at R = 1.0 for W are introduced. In proof of this, 
a case was run which took’as an initial state the infinite disk solution. 


Convergence still required 2000 iterations. 


In order to test the stability of the numerical solution, two 
special computer runs were made. “In the first, a convergent solution 
was perturbed in W. This perturbation was of the order ae |w| and 
was introduced between R = 0.5 and 0.6. The solution ‘gain con- 
verged, suggesting stability for finite perturbations. For the second 
run, initial values were assumed from the infinite disk solution and 


pressure was increased by an order of magnitude. This solution 


achieved an oscillatory behavior then diverged completely. 
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VI. CONCLUSIONS 


A large scale parametric examination of §, AT, AR and AZ is 
required to best determine the numerical stability, the rate of con- 
vergence, and the accuracy of each problem solved using Chorin's 
technique. To show that the venules obtained herein are indeed 
reasonable a numerical solution is compared with a series solution 
for infinite disks. The series solutions are written to second order 


in 2X R, as follows [ 2, 3]: 


U=- RAR, [4 (1-Z) - 9 (1-2)? +5 (1-z)* J 
to 


V = R (1-2) - RA*R.” (-8 (1-Z) - 35 (1-Z)* + 63 (1-Z)®) - 20 (1-Z)’ J 
6300 


¢ 


W=AR, (- 4 (1-2)? + 6 (1-2)? - 2 (1-Z)® J 
60 
P=.15R®, 


These comparisons are made for \ = .1 and R= 5.0. Figure 7 shows 


the secondary velocity field for the numerical solution. (A computer 


output for this solution is contained in Appendix B. ) 
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T 


Scale: .25" = . 00369 


ref 


R=1.0 


Figure 8 shows a graphical representation of a comparison between 
the numerical solution for finite disks and the serics solution for 


infinite disks. 


Considering the closeness of the results of the above compari- 
son it can definitely be stated that the velocity and pressure for 
rotating diskviscometer applications are adequately represented by 


the series solution of the infinite disk problem. 
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APPENDIX A 


COMPUTER PROGRAM ' 


Coding hints: 


1. Factor all constants and remove their computation from 
the innermost integration loop. This results in a 30% 


+ 


reduction of computer time. © 


2. Simplify the subscript arithmetic as much as possible, 
(See statement number 5.) This results ina 20% reduction 


of computer time. 
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APPENDIX B 


SAMPLE COMPUTER OUTPUTS 


Guide to the computer output: 
1. The first line of print is the number of iterations. 


2. The second line of print is A," Ro 6, AT, AR, AZ anda 


. pair of indices. 


¢ 


3, The'velocity fields (U, V, W) and the density (D=6P) are 
printed on the following lines. The numbers are oriented 


so as to correspond to the graphs previously shown. 
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